/********************************************************************
File contains routine:
  apr_gauss_select - to select the proper quadrature
  apr_gauss_init - to initialize Gauss-Legendre quadratures

------------------------------
History:
	05.2001 - Krzysztof Banas, initial version
*********************************************************************/

#include<stdlib.h>
#include<stdio.h>
#include<math.h>

/* internal header file for the standard approximation module */
#include "./aph_std_lin.h"

static double Xg_1d[4][4];
static double Wg_1d[4][4];
static double Xg_2d[5][21];
static double Wg_2d[5][7];
static double Xg_3d[5][60];
static double Wg_3d[5][15];



/*--------------------------------------------------------------------
  apr_gauss_select - to select the proper quadrature
--------------------------------------------------------------------*/
void apr_gauss_select(
	int Type,	/* in: type of integration considered */
			/* 	1 - LINE INTEGRATION over [-1,1] */
			/* 	2 - GAUSSIAN QUADRATURES FOR TRIANGLES */
	                /*      3 - GAUSSAIN QUADRATURES FOR TETRAHEDRONS */
	int Order, 	/* in: order of integration == degree */
			/*     of polynomials integrated exactly */
        int *Ng,	/* out: number of gauss points */
	double **Xg_p, 	/* out: gauss points */
	double **Wg_p 	/* out: gauss weights */
	)
{
	/* auxiliary variables */
//	int i,j;

	/*++++++++++++++++ executable statements ++++++++++++++++*/

  switch(Type)
    {
    case 1: 
      { /* LINE INTEGRATION over [-1,1] */
	*Ng=Order/2+1;
	
	if(*Ng<=4){
	  *Xg_p = Xg_1d[*Ng - 1];
	  *Wg_p = Wg_1d[*Ng - 1];
	}
	else {
	  printf("not implemented yet, number of gauss points too big!\n");
	  exit(-1);
	}
      } 
      break;
    case 2:
      { /* GAUSSIAN QUADRATURES FOR TRIANGLES */
	switch(Order)
	  {
	  case 0: case 1: *Ng  = 1; break;
	  case 2: case 3: *Ng = Order+1; break;
	  case 4: case 5: *Ng = Order+2; break;
	  }
	
	if(Order<=0){
	  *Xg_p = Xg_2d[0];
	  *Wg_p = Wg_2d[0];
	}
	else if(Order<=5){
	  *Xg_p = Xg_2d[Order - 1];
	  *Wg_p = Wg_2d[Order - 1];
	  
	}
	else {
	  printf("not implemented yet, number of gauss points too big!\n");
	  exit(-1);
	}
      }
      break;
    case 3: 
      { //  3 - GAUSSAIN QUADRATURES FOR TETRAHEDRONS
	switch(Order)
	  {
	  case 0: case 1: *Ng  = 1; break;
	  case 2: case 3: *Ng = Order+2; break;
	  case 4: *Ng = 11; break;
	  case 5: *Ng = 15; break;
	  }
	
	if(Order<=0){
	  *Xg_p = Xg_3d[0];
	  *Wg_p = Wg_3d[0];
	}
	else if(Order<=5){
	  *Xg_p = Xg_3d[Order - 1];
	  *Wg_p = Wg_3d[Order - 1];
	}
      }
      break;
    default:
      {
	printf("not implemented yet, number of gauss points too big!\n");
	exit(-1);
      } 
      break;
    }
  
}

/*--------------------------------------------------------------------
  apr_gauss_init - to initialize Gauss-Legendre quadratures
--------------------------------------------------------------------*/
void apr_gauss_init()
{

  /* auxiliary variables */
  int i,j;

/*++++++++++++++++ executable statements ++++++++++++++++*/

  for(i=0;i<4;i++){
    for(j=0;j<4;j++){
      Xg_1d[i][j]=0.0;
      Wg_1d[i][j]=0.0;
    }
  }

  for(i=0;i<5;i++){
    for(j=0;j<21;j++){
      Xg_2d[i][j]=0.0;
    }
    for(j=0;j<7;j++){
      Wg_2d[i][j]=0.0;
    }
  }

  for(i=0;i<5;i++){
	for(j=0;j<60;j++){
	  Xg_3d[i][j]=0.0;
	}
	for(j=0;j<15;j++){
	  Wg_3d[i][j]=0.0;
	}
  }

	Xg_1d[0][0] = 0.0;
	Wg_1d[0][0] = 2.0;

	Xg_1d[1][0] = -0.577350269189626;
	Xg_1d[1][1] =  0.577350269189626;
	Wg_1d[1][0] = 1.0;
	Wg_1d[1][1] = 1.0;

	Xg_1d[2][0] = -0.774596669241483;
	Xg_1d[2][1] =  0.000000000000000;
	Xg_1d[2][2] =  0.774596669241483;
	Wg_1d[2][0] =  0.555555555555553;
	Wg_1d[2][1] =  0.888888888888889;
	Wg_1d[2][2] =  0.555555555555553;

	Xg_1d[3][0] = -0.861136311594053;
	Xg_1d[3][1] = -0.339981043584856;
	Xg_1d[3][2] =  0.339981043584856;
	Xg_1d[3][3] =  0.861136311594053;
	Wg_1d[3][0] =  0.347854845137448;
	Wg_1d[3][1] =  0.652145154862546;
	Wg_1d[3][2] =  0.652145154862546;
	Wg_1d[3][3] =  0.347854845137448;


	Wg_2d[0][ 0]=   1.000000000000000;
	Xg_2d[0][  0]=   0.333333333333333;
	Xg_2d[0][  1]=   0.333333333333333;
	Xg_2d[0][  2]=   0.333333333333333;

	Wg_2d[1][ 0]=   0.333333333333333;
	Xg_2d[1][  0]=   0.666666666666667;
	Xg_2d[1][  1]=   0.166666666666667;
	Xg_2d[1][  2]=   0.166666666666667;
	Wg_2d[1][ 1]=   0.333333333333333;
	Xg_2d[1][  3]=   0.166666666666667;
	Xg_2d[1][  4]=   0.166666666666667;
	Xg_2d[1][  5]=   0.666666666666667;
	Wg_2d[1][ 2]=   0.333333333333333;
	Xg_2d[1][  6]=   0.166666666666667;
	Xg_2d[1][  7]=   0.666666666666667;
	Xg_2d[1][  8]=   0.166666666666667;

	Wg_2d[2][ 0]=  -0.562500000000000;
	Xg_2d[2][  0]=   0.333333333333333;
	Xg_2d[2][  1]=   0.333333333333333;
	Xg_2d[2][  2]=   0.333333333333333;
	Wg_2d[2][ 1]=   0.520833333333333;
	Xg_2d[2][  3]=   0.600000000000000;
	Xg_2d[2][  4]=   0.200000000000000;
	Xg_2d[2][  5]=   0.200000000000000;
	Wg_2d[2][ 2]=   0.520833333333333;
	Xg_2d[2][  6]=   0.200000000000000;
	Xg_2d[2][  7]=   0.200000000000000;
	Xg_2d[2][  8]=   0.600000000000000;
	Wg_2d[2][ 3]=   0.520833333333333;
	Xg_2d[2][  9]=   0.200000000000000;
	Xg_2d[2][ 10]=   0.600000000000000;
	Xg_2d[2][ 11]=   0.200000000000000;

	Wg_2d[3][ 0]=   0.223381589678011;
	Xg_2d[3][  0]=   0.108103018168070;
	Xg_2d[3][  1]=   0.445948490915965;
	Xg_2d[3][  2]=   0.445948490915965;
	Wg_2d[3][ 1]=   0.223381589678011;
	Xg_2d[3][  3]=   0.445948490915965;
	Xg_2d[3][  4]=   0.445948490915965;
	Xg_2d[3][  5]=   0.108103018168070;
	Wg_2d[3][ 2]=   0.223381589678011;
	Xg_2d[3][  6]=   0.445948490915965;
	Xg_2d[3][  7]=   0.108103018168070;
	Xg_2d[3][  8]=   0.445948490915965;
	Wg_2d[3][ 3]=   0.109951743655322;
	Xg_2d[3][  9]=   0.816847572980459;
	Xg_2d[3][ 10]=   0.091576213509771;
	Xg_2d[3][ 11]=   0.091576213509771;
	Wg_2d[3][ 4]=   0.109951743655322;
	Xg_2d[3][ 12]=   0.091576213509771;
	Xg_2d[3][ 13]=   0.091576213509771;
	Xg_2d[3][ 14]=   0.816847572980459;
	Wg_2d[3][ 5]=   0.109951743655322;
	Xg_2d[3][ 15]=   0.091576213509771;
	Xg_2d[3][ 16]=   0.816847572980459;
	Xg_2d[3][ 17]=   0.091576213509771;

	Wg_2d[4][ 0]=   0.225000000000000;
	Xg_2d[4][  0]=   0.333333333333333;
	Xg_2d[4][  1]=   0.333333333333333;
	Xg_2d[4][  2]=   0.333333333333333;
	Wg_2d[4][ 1]=   0.132394152788506;
	Xg_2d[4][  3]=   0.059715871789770;
	Xg_2d[4][  4]=   0.470142064105115;
	Xg_2d[4][  5]=   0.470142064105115;
	Wg_2d[4][ 2]=   0.132394152788506;
	Xg_2d[4][  6]=   0.470142064105115;
	Xg_2d[4][  7]=   0.470142064105115;
	Xg_2d[4][  8]=   0.059715871789770;
	Wg_2d[4][ 3]=   0.132394152788506;
	Xg_2d[4][  9]=   0.470142064105115;
	Xg_2d[4][ 10]=   0.059715871789770;
	Xg_2d[4][ 11]=   0.470142064105115;
	Wg_2d[4][ 4]=   0.125939180544827;
	Xg_2d[4][ 12]=   0.797426985353087;
	Xg_2d[4][ 13]=   0.101286507323456;
	Xg_2d[4][ 14]=   0.101286507323456;
	Wg_2d[4][ 5]=   0.125939180544827;
	Xg_2d[4][ 15]=   0.101286507323456;
	Xg_2d[4][ 16]=   0.101286507323456;
	Xg_2d[4][ 17]=   0.797426985353087;
	Wg_2d[4][ 6]=   0.125939180544827;
	Xg_2d[4][ 18]=   0.101286507323456;
	Xg_2d[4][ 19]=   0.797426985353087;
	Xg_2d[4][ 20]=   0.101286507323456;
	// p=1,n=1
	Wg_3d[0][0]=1.0;
	Xg_3d[0][ 0]=0.25;
	Xg_3d[0][ 1]=0.25;
	Xg_3d[0][ 2]=0.25;
	Xg_3d[0][ 3]=0.25;
	// p=2,n=4
	Wg_3d[1][0]=0.25;
	Xg_3d[1][ 0]=0.585410196624969;
	Xg_3d[1][ 1]=0.138196601125011;
	Xg_3d[1][ 2]=0.138196601125011;
	Xg_3d[1][ 3]=0.138196601125011;
	Wg_3d[1][1]=0.25;
	Xg_3d[1][ 4]=0.138196601125011;
	Xg_3d[1][ 5]=0.585410196624969;
	Xg_3d[1][ 6]=0.138196601125011;
	Xg_3d[1][ 7]=0.138196601125011;
	Wg_3d[1][2]=0.25;
	Xg_3d[1][ 8]=0.138196601125011;
	Xg_3d[1][ 9]=0.138196601125011;
	Xg_3d[1][ 10]=0.585410196624969;
	Xg_3d[1][ 11]=0.138196601125011;
	Wg_3d[1][3]=0.25;
	Xg_3d[1][ 12]=0.138196601125011;
	Xg_3d[1][ 13]=0.138196601125011;
	Xg_3d[1][ 14]=0.138196601125011;
	Xg_3d[1][ 15]=0.585410196624969;
	//p=3,n=5
	Wg_3d[2][0]=-0.8;
	Xg_3d[2][ 0]=0.25;
	Xg_3d[2][ 1]=0.25;
	Xg_3d[2][ 2]=0.25;
	Xg_3d[2][ 3]=0.25;
	Wg_3d[2][1]=0.45;
	Xg_3d[2][ 4]=0.5;
	Xg_3d[2][ 5]=0.1666666666666667;
	Xg_3d[2][ 6]=0.1666666666666667;
	Xg_3d[2][ 7]=0.1666666666666667;
	Wg_3d[2][2]=0.45;
	Xg_3d[2][ 8]=0.1666666666666667;
	Xg_3d[2][ 9]=0.5;
	Xg_3d[2][10]=0.1666666666666667;
	Xg_3d[2][11]=0.1666666666666667;
	Wg_3d[2][3]=0.45;
	Xg_3d[2][ 8]=0.1666666666666667;
	Xg_3d[2][ 9]=0.1666666666666667;
	Xg_3d[2][10]=0.5;
	Xg_3d[2][11]=0.1666666666666667;
	Wg_3d[2][4]=0.45;
	Xg_3d[2][12]=0.1666666666666667;
	Xg_3d[2][13]=0.1666666666666667;
	Xg_3d[2][14]=0.1666666666666667;
	Xg_3d[2][15]=0.5;
	// p=4,n=11
	Wg_3d[3][0]=-0.013155555555556;
	Xg_3d[3][ 0]=0.25;
	Xg_3d[3][ 1]=0.25;
	Xg_3d[3][ 2]=0.25;
	Xg_3d[3][ 3]=0.25;
	Wg_3d[3][1]=0.007622222222222;
	Xg_3d[3][ 4]=0.785714285714286;
	Xg_3d[3][ 5]=0.071428571428571;
	Xg_3d[3][ 6]=0.071428571428571;
	Xg_3d[3][ 7]=0.071428571428571;
	Wg_3d[3][2]=0.007622222222222;
	Xg_3d[3][ 8]=0.071428571428571;
	Xg_3d[3][ 9]=0.785714285714286;
	Xg_3d[3][10]=0.071428571428571;
	Xg_3d[3][11]=0.071428571428571;
	Wg_3d[3][3]=0.007622222222222;
	Xg_3d[3][12]=0.071428571428571;
	Xg_3d[3][13]=0.071428571428571;
	Xg_3d[3][14]=0.785714285714286;
	Xg_3d[3][15]=0.071428571428571;
	Wg_3d[3][4]=0.007622222222222;
	Xg_3d[3][16]=0.071428571428571;
	Xg_3d[3][17]=0.071428571428571;
	Xg_3d[3][18]=0.071428571428571;
	Xg_3d[3][19]=0.785714285714286;
	Wg_3d[3][5]=0.024888888888889;
	Xg_3d[3][20]=0.399403576166799;
	Xg_3d[3][21]=0.100596423833201;
	Xg_3d[3][22]=0.399403576166799;
	Xg_3d[3][23]=0.100596423833201;
	Wg_3d[3][6]=0.024888888888889;
	Xg_3d[3][24]=0.100596423833201;
	Xg_3d[3][25]=0.399403576166799;
	Xg_3d[3][26]=0.100596423833201;
	Xg_3d[3][27]=0.399403576166799;
	Wg_3d[3][7]=0.024888888888889;
	Xg_3d[3][28]=0.399403576166799;
	Xg_3d[3][29]=0.399403576166799;
	Xg_3d[3][30]=0.100596423833201;
	Xg_3d[3][31]=0.100596423833201;
	Wg_3d[3][8]=0.024888888888889;
	Xg_3d[3][32]=0.100596423833201;
	Xg_3d[3][33]=0.100596423833201;
	Xg_3d[3][34]=0.399403576166799;
	Xg_3d[3][35]=0.399403576166799;
	Wg_3d[3][9]=0.024888888888889;
	Xg_3d[3][36]=0.100596423833201;
	Xg_3d[3][37]=0.399403576166799;
	Xg_3d[3][38]=0.399403576166799;
	Xg_3d[3][39]=0.100596423833201;
	Wg_3d[3][10]=0.024888888888889;
	Xg_3d[3][40]=0.399403576166799;
	Xg_3d[3][41]=0.100596423833201;
	Xg_3d[3][42]=0.100596423833201;
	Xg_3d[3][43]=0.399403576166799;
	// p=5, n=15
	Wg_3d[4][ 0]=0.030283678097089;
	Xg_3d[4][0]=0.25;
	Xg_3d[4][1]=0.25;
	Xg_3d[4][2]=0.25;
	Xg_3d[4][3]=0.25;

	Wg_3d[4][ 1]=0.006026785714286;
	Xg_3d[4][4]=0.0;
	Xg_3d[4][5]=0.333333333333333;
	Xg_3d[4][6]=0.333333333333333;
	Xg_3d[4][7]=0.333333333333333;
	Wg_3d[4][2]=0.006026785714286;
	Xg_3d[4][8]=0.333333333333333;
	Xg_3d[4][ 9]=0.0;
	Xg_3d[4][10]=0.333333333333333;
	Xg_3d[4][11]=0.333333333333333;
	Wg_3d[4][ 3]=0.006026785714286;
	Xg_3d[4][12]=0.333333333333333;
	Xg_3d[4][13]=0.333333333333333;
	Xg_3d[4][14]=0.0;
	Xg_3d[4][15]=0.333333333333333;
	Wg_3d[4][ 4]=0.006026785714286;
	Xg_3d[4][16]=0.333333333333333;
	Xg_3d[4][17]=0.333333333333333;
	Xg_3d[4][18]=0.333333333333333;
	Xg_3d[4][19]=0.0;

	Wg_3d[4][ 5]=0.011645249086029;
	Xg_3d[4][20]=0.727272727272727;
	Xg_3d[4][21]=0.090909090909091;
	Xg_3d[4][22]=0.090909090909091;
	Xg_3d[4][23]=0.090909090909091;
	Wg_3d[4][ 6]=0.011645249086029;
	Xg_3d[4][24]=0.090909090909091;
	Xg_3d[4][25]=0.727272727272727;
	Xg_3d[4][26]=0.090909090909091;
	Xg_3d[4][27]=0.090909090909091;
	Wg_3d[4][ 7]=0.011645249086029;
	Xg_3d[4][28]=0.090909090909091;
	Xg_3d[4][29]=0.090909090909091;
	Xg_3d[4][30]=0.727272727272727;
	Xg_3d[4][31]=0.090909090909091;
	Wg_3d[4][ 8]=0.011645249086029;
	Xg_3d[4][32]=0.090909090909091;
	Xg_3d[4][33]=0.090909090909091;
	Xg_3d[4][34]=0.090909090909091;
	Xg_3d[4][35]=0.727272727272727;

	Wg_3d[4][ 9]=0.010949141561386 ;
	Xg_3d[4][36]=0.066550153573664;
	Xg_3d[4][37]=0.433449846426336;
	Xg_3d[4][38]=0.066550153573664;
	Xg_3d[4][39]=0.433449846426336;
		Wg_3d[4][10]=0.010949141561386 ;
	Xg_3d[4][40]=0.433449846426336;
	Xg_3d[4][41]=0.066550153573664;
	Xg_3d[4][42]=0.433449846426336;
	Xg_3d[4][43]=0.066550153573664;
		Wg_3d[4][11]=0.010949141561386 ;
	Xg_3d[4][44]=0.066550153573664;
	Xg_3d[4][45]=0.066550153573664;
	Xg_3d[4][46]=0.433449846426336;
	Xg_3d[4][47]=0.433449846426336;
		Wg_3d[4][12]=0.010949141561386 ;
	Xg_3d[4][48]=0.433449846426336;
	Xg_3d[4][49]=0.433449846426336;
	Xg_3d[4][50]=0.066550153573664;
	Xg_3d[4][51]=0.066550153573664;
		Wg_3d[4][13]=0.010949141561386 ;
	Xg_3d[4][52]=0.066550153573664;
	Xg_3d[4][53]=0.433449846426336;
	Xg_3d[4][54]=0.433449846426336;
	Xg_3d[4][55]=0.066550153573664;
		Wg_3d[4][14]=0.010949141561386 ;
	Xg_3d[4][56]=0.433449846426336;
	Xg_3d[4][57]=0.066550153573664;
	Xg_3d[4][58]=0.066550153573664;
	Xg_3d[4][59]=0.433449846426336;

}

